{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "ec4d5f74",
   "metadata": {},
   "source": [
    "# 机器读心术之文本挖掘与自然语言处理第4课书面作业\n",
    "\n",
    "学号：207402  \n",
    "\n",
    "**书面作业**\n",
    "\n",
    "1. 《概率论与数理统计》第333页第7题  \n",
    "2. 体验HMM第一个基本问题前向算法或后向算法解法，计算本课程幻灯片第34页“掷骰子HMM”中出现可观测序列为“1635273524”的概率，可使用任何编程语言或手工计算"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "318c00a6",
   "metadata": {},
   "source": [
    "## 作业1\n",
    "设任意相继的两天中，雨天转晴天的概率为$\\frac{1}{3}$，晴天转雨天的概率为$\\frac{1}{2}$，任一天晴或雨是互为逆事件。以0表示晴天状态，以1表示雨天状态，$X_n$表示第n天的状态（0或1）。试写出马氏链$\\{X_n, n \\ge 1\\}$的一步转移概率矩阵。又若已经5月1日为晴天，问5月3日为晴天，5月5日为雨天的概率各是多少？"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "78c94acb",
   "metadata": {},
   "source": [
    "**解：**  \n",
    "一步转换概率矩阵如下：  \n",
    "第一行为雨天，第二行为晴天。  \n",
    "第一列为雨天，第二列为晴天。"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "90e7e150",
   "metadata": {},
   "source": [
    "$$\n",
    "T=\n",
    "\\left[\n",
    "\\begin{matrix}\n",
    "\\frac{2}{3} & \\frac{1}{3}\\\\\n",
    "\\frac{1}{2} & \\frac{1}{2} \\\\\n",
    "\\end{matrix}\n",
    "\\right]\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "c41db141",
   "metadata": {},
   "source": [
    "$X_0$表示5月1日状态，$X_3$表示5月3日状态，$X_5$表示5月5日状态，则："
   ]
  },
  {
   "cell_type": "markdown",
   "id": "8713b64f",
   "metadata": {},
   "source": [
    "$$\n",
    "X_0=\n",
    "\\left[\n",
    "\\begin{matrix}\n",
    "0 & 1 \\\\\n",
    "\\end{matrix}\n",
    "\\right] \\\\\n",
    "X_3=X_0 * T * T \\\\\n",
    "X_5=X_3 * T * T\n",
    "$$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "id": "b32acb8f",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "X3= [0.58333333 0.41666667]\n",
      "X5= [0.59953704 0.40046296]\n",
      "5月3日为晴天的概率： 0.41666666666666663\n",
      "5月5日为雨天的概率： 0.5995370370370369\n"
     ]
    }
   ],
   "source": [
    "import numpy as np\n",
    "\n",
    "T=np.array(\n",
    "    [\n",
    "        [2/3.,1/3.],\n",
    "        [1/2.,1/2.]]\n",
    ")\n",
    "\n",
    "X0=np.array([0,1])\n",
    "X3=np.matmul(np.matmul(X0,T),T)\n",
    "X5=np.matmul(np.matmul(X3,T),T)\n",
    "print('X3=',X3)\n",
    "print('X5=',X5)\n",
    "print('5月3日为晴天的概率：',X3[1])\n",
    "print('5月5日为雨天的概率：',X5[0])"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "f4e494fe",
   "metadata": {},
   "source": [
    "## 作业2"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "955b75a1",
   "metadata": {},
   "source": [
    "假设手里有三个不同的骰子。第一个骰子是我们平常见的骰子（称这个骰子为D6）， 6个面，每个面（1，2，3，4，5，6）出现的概率是1/6。第二个骰子是个四面体（称这个骰子为D4），每个面（1，2，3，4）出现的概率是1/4。第三个骰子有八个面 （称这个骰子为D8），每个面（1，2，3，4，5，6，7，8）出现的概率是1/8。"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "9cb290c8",
   "metadata": {},
   "source": [
    "![nlp04-1](https://i.loli.net/2021/11/09/WOI7nSyxPEcjp8s.png)\n",
    "\n",
    "![nlp04-2](https://i.loli.net/2021/11/09/ah8BbqCOoWKPIGF.png)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "ece804bf",
   "metadata": {},
   "source": [
    "**解：**  \n",
    "编码如下。"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "id": "23dbf70d",
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "\n",
    "# 状态转换概率矩阵定义如下：\n",
    "a=np.array(\n",
    "    [\n",
    "        [1/3.,1/3.,1/3.],\n",
    "        [1/3.,1/3.,1/3.],\n",
    "        [1/3.,1/3.,1/3.]\n",
    "    ]\n",
    ")\n",
    "\n",
    "# 发射概率矩阵定义如下：\n",
    "b=np.array(\n",
    "    [\n",
    "        [1/6,1/6,1/6,1/6,1/6,1/6,  0,  0],\n",
    "        [1/4,1/4,1/4,1/4,  0,  0,  0,  0],\n",
    "        [1/8,1/8,1/8,1/8,1/8,1/8,1/8,1/8]\n",
    "    ]\n",
    ")\n",
    "\n",
    "#初始概率分布定义如下：\n",
    "pi=np.array([1/3,1/3,1/3])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "id": "4079a165",
   "metadata": {},
   "outputs": [],
   "source": [
    "def forward_proc(a,b,pi,seq,N):\n",
    "    '''\n",
    "    前向算法\n",
    "    a: 状态转换概率矩阵\n",
    "    b: 发射概率矩阵\n",
    "    pi: 初始概率分布\n",
    "    seq: 观察序列\n",
    "    N: 状态数\n",
    "    '''\n",
    "    T=len(seq)\n",
    "    alpha=[]\n",
    "    alpha.append(pi*b[:,int(seq[0])])\n",
    "    for t in range(T-1):\n",
    "        alpha_t=alpha[-1]\n",
    "        alpha_t_plus_1=np.matmul(alpha_t,a)*b[:,int(seq[t+1])]\n",
    "        alpha.append(alpha_t_plus_1)\n",
    "    return np.sum(alpha[-1])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "id": "e61fdd15",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "3.0614588526568534e-10\n"
     ]
    }
   ],
   "source": [
    "print(forward_proc(a,b,pi,'1635273524',3))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "id": "72eff8d7",
   "metadata": {},
   "outputs": [],
   "source": [
    "def backward_proc(a,b,pi,seq,N):\n",
    "    '''\n",
    "    后向算法\n",
    "    a: 状态转换概率矩阵\n",
    "    b: 发射概率矩阵\n",
    "    pi: 初始概率分布\n",
    "    seq: 观察序列\n",
    "    N: 状态数\n",
    "    '''\n",
    "    T=len(seq)\n",
    "    beta=[]\n",
    "    beta.append([1. for i in range(N)])\n",
    "    for t in range(T-1,0,-1):\n",
    "        beta_t_plus_1=beta[-1]\n",
    "        beta_t_plus_1=np.matmul(a,b[:,int(seq[t+1-1])])*beta_t_plus_1\n",
    "        beta.append(beta_t_plus_1)\n",
    "    return np.sum(pi*b[:,int(seq[0])]*beta[-1])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "id": "c54f934b",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "3.0614588526568513e-10\n"
     ]
    }
   ],
   "source": [
    "print(backward_proc(a,b,pi,'1635273524',3))"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.8.8"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}
